package scu.maqiang.electromagnetics;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.*;

import java.util.Arrays;

/***
 * @author 马强
 * @brief 平行板电容器问题
 *        计算平行板电容器板内与周围的电势与电场分部
 *        例子来源于MATLAB®-based Finite Element Programming in Electromagnetic Modeling
 */
public class ParallelPlateCapacitor {
    public static void main(String[] args) {
        //parallelPlateWithTrancatedComputationalDomain();
        //parallelPlateWithDielectricInhomogeneity();
        microstripTransmissionLine();
    }

    public static void parallelPlateWithTrancatedComputationalDomain() {
        String meshFileName = "C:\\Users\\Administrator\\OneDrive\\FFStudy\\Homogenization\\electromagnetics\\ParallelPlate.msh";
        Mesh2T3 mesh = new Mesh2T3();
        mesh.readMesh(MeshFormat.msh, meshFileName);
        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double eps0 = 1.0;
        double epsd = 1.0 * eps0;
        ScalarFunc coef = (xy, label, param) -> label == 4? eps0: epsd;
        fs.assembleStiff(new ScalarFunc[]{coef}, null, BVPType.COMMON, A);
        fs.applyBC_MBN(A, RHS, Direct.All, 0.0, 1, 2, 3, 4, 11);
        fs.applyBC_MBN(A, RHS, Direct.All, 1.0, 13);
        double[][] grad = new double[2][fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver(A);
        double[] x = new double[fs.GetNdof()];
        solver.PCGSSOR(RHS, x, 1.5, 1);
        fs.computeGradient(x, grad);
        Arrays.setAll(grad[0], i->-grad[0][i]);
        Arrays.setAll(grad[1], i->-grad[1][i]);
        System.out.println("Solution:");
        System.out.println("min: " + MVO.min(x) + "\t\t" + "max: " + MVO.max(x));
        mesh.toTecplot("Potential_ParallelPlate.dat", x);
        mesh.toTecplot("EField_ParallelPlate.dat", grad);
        
    }

    public static void parallelPlateWithDielectricInhomogeneity() {
        String meshFileName = "C:\\Users\\Administrator\\OneDrive\\FFStudy\\Homogenization\\electromagnetics\\ParallelPlate2.msh";
        Mesh2T3 mesh = new Mesh2T3();
        mesh.readMesh(MeshFormat.msh, meshFileName);
        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double eps0 = 1.0;
        double epsd = 16.0 * eps0;
        ScalarFunc coef = (xy, label, param) -> label == 4? eps0: epsd;
        fs.assembleStiff(new ScalarFunc[]{coef}, null, BVPType.COMMON, A);
        //A.sort();
        System.out.println(A.getNNZ());
        //System.out.println(A.display());
        //.applyBC_RCE(A, RHS, Direct.All, 1.0, 3);
        //fs.applyBC_RCE(A, RHS, Direct.All, -1.0, 1);
        fs.applyBC_MBN(A, RHS, Direct.All, 1.0, 3);
        fs.applyBC_MBN(A, RHS, Direct.All, -1.0, 1);
        double[][] grad = new double[2][fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver(A);
        double[] x = new double[fs.GetNdof()];
        //IterSSolver solver = new IterSSolver();
        //solver.PCGSSOR(RHS, x, 1.5, 1);
        solver.CG(RHS, x, 2);
        //solver.PCGSSOR(A, RHS, x,1.5, 1);
        fs.computeGradient(x, grad);
        Arrays.setAll(grad[0], i->-grad[0][i]);
        Arrays.setAll(grad[1], i->-grad[1][i]);
        System.out.println("Solution:");
        System.out.println("min: " + MVO.min(x) + "\t\t" + "max: " + MVO.max(x));
        mesh.toTecplot("Potential_ParallelPlate2.dat", x);
        mesh.toTecplot("EField_ParallelPlate2.dat", grad);
    }

    public static void microstripTransmissionLine() {
        String meshFileName = "C:\\Users\\Administrator\\OneDrive\\FFStudy\\Homogenization\\electromagnetics\\microstripLine.msh";
        Mesh2T3 mesh = new Mesh2T3();
        mesh.readMesh(MeshFormat.msh, meshFileName);
        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double eps0 = 8.85e-12;
        double epsd = 4.0 * eps0;
        ScalarFunc coef = (xy, label, param) -> label == 4? eps0: epsd;
        fs.assembleStiff(new ScalarFunc[]{coef}, null, BVPType.COMMON, A);
        //A.sort();
        //System.out.println(A.getNNZ());
        //System.out.println(A.display());
        //.applyBC_RCE(A, RHS, Direct.All, 1.0, 3);
        //fs.applyBC_RCE(A, RHS, Direct.All, -1.0, 1);
        fs.applyBC_MBN(A, RHS, Direct.All, 0.0, 1, 21, 22, 3, 41, 42);
        fs.applyBC_MBN(A, RHS, Direct.All, 1.0, 52);
        double[][] grad = new double[2][fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver(A);
        double[] x = new double[fs.GetNdof()];
        //IterSSolver solver = new IterSSolver();
        solver.PCGSSSOR(RHS, x, 1.5, 1);
        //solver.CG(RHS, x, 2);
        //solver.PCGSSOR(A, RHS, x,1.5, 1);
        fs.computeGradient(x, grad);
        Arrays.setAll(grad[0], i->-grad[0][i]);
        Arrays.setAll(grad[1], i->-grad[1][i]);
        System.out.println("Solution:");
        System.out.println("min: " + MVO.min(x) + "\t\t" + "max: " + MVO.max(x));
        mesh.toTecplot("Potential_microstripLine.dat", x);
        mesh.toTecplot("EField_microstripLine.dat", grad);

    }


}
